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A Parallel Dual Fast Gradient Method for MPC Applications* 


Laura Ferranti^, Tamas Keviczky^ 


Abstract —We propose a parallel adaptive constraint¬ 
tightening approach to solve a linear model predictive control 
problem for discrete-time systems, based on inexact numerical 
optimization algorithms and operator splitting methods. The 
underlying algorithm first splits the original problem in as 
many independent subproblems as the length of the prediction 
horizon. Then, our algorithm computes a solution for these 
subproblems in parallel by exploiting auxiliary tightened sub¬ 
problems in order to certify the control law in terms of subopti¬ 
mality and recursive feasibility, along with closed-loop stability 
of the controlled system. Compared to prior approaches based 
on constraint tightening, our algorithm computes the tightening 
parameter for each subproblem to handle the propagation 
of errors introduced by the parallelization of the original 
problem. Our simulations show the computational benefits 
of the parallelization with positive impacts on performance 
and numerical conditioning when compared with a recent 
nonparallel adaptive tightening scheme. 

I. Introduction 

Model Predictive Control (MPC) is a consolidated control 
technique that can efficiently handle constraints on the pro¬ 
cess to be controlled. Nevertheless, its application is not yet 
widespread in many domains where real-time computational 
constraints and requirements of certihed solutions are of 
major concern, such as aerospace or automotive applications. 
There is a growing interest in both industry and academia for 
exploring parallel solutions to MPC problems ([1], [2], [3]), 
especially in light of the emerging many-core architectures, 
aiming to improve the computational efficiency of solving 
the underlying optimization problem. 

Contribution. In this paper, we explore the use of pa¬ 
rallelization techniques to efficiently solve a typical MPC 
problem for a linear discrete-time system, with a substan¬ 
tial computational speedup compared to nonparallel imple¬ 
mentations. Our proposed algorithm combines the use of 
Alternating Direction Method of Multipliers (ADMMs [4], 
[5]) to handle the coupling constraints that arise from the 
dynamics of the system and inexact solvers (i.e., solvers 
that guarantee feasibility and optimality only asymptotically 
with the number of iterations), such as the Nesterov’s Dual 
Fast Gradient (DFG) method [10]. In particular, the first 
step of the proposed algorithm is to split the original MPC 
problem over the length N of the prediction horizon into 
TV 4-1 independent subproblems {time-splitting [3]) solved by 
N -\-l parallel workers periodically exchanging information 
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at predetermined synchronization points. Then, the second 
step is to solve these subproblems in parallel using an 
inexact solver and guarantee, at the same time, that the 
solution of the original MPC problem is recursively feasible 
and the system is closed-loop stable. The combination of 
parallelization and inexact solvers can result in infeasibility 
and closed-loop instability. We rely on an algorithm based 
on constraint tightening to overcome these issues. Loosely 
speaking, constraint-tightening algorithms solve an alterna¬ 
tive problem in which the constraints have been tightened 
by a certain amount to compensate for the accuracy loss 
(and possible related infeasibility) introduced by the solver. 
We rely on an adaptive tightening strategy to select an 
appropriate amount of tightening for our algorithm. Every 
time new measurements are available from the plant, our 
algorithm chooses the amount of tightening required for each 
subproblem in order to compensate for the error introduced 
by the time-splitting combined with the inexact solver. 

Related work. The time-splitting technique has been 
proposed in [3]. In contrast to [3], we combine ADMM with 
inexact solvers and focus on the requirements for recursive 
feasibility and closed-loop stability of the original problem. 

Other constraint-tightening schemes have been proposed in 
the literature (outside the parallel framework). For example, 
the authors in [7] propose an algorithm in which the amount 
of tightening is chosen offline to guarantee suboptimality 
and feasibility of the solution for all the initial states of 
the MPC problem. Solutions based on adaptive constraint 
tightening have been recently proposed in [8], where the 
tightening parameter is chosen adaptively. Compared to [8], 
our tightening update rule allows for a nonuniform amount 
of tightening (the tightening varies along the prediction 
horizon). Furthermore, thanks to the modular structure of 
our approach, the optimizer solves simpler problems of hxed 
dimension, which is independent from N. As a consequence, 
an increase of N does not affect the conditioning of the 
problem and the convergence of the solver. Hence, our 
approach leads to a performance improvement even when 
forcing full serialization of the parallel operations (i.e., 
serialized mode [9]). 

Outline. In the following. Section |I^ presents the initial 
problem formulation. Section III introduces the auxiliary 
subproblems and our proposed solver. Section |IV] describes 
our strategy to select the tightening of each subproblem to 
handle the parallelization error. Section [V]proposes an online 
update strategy of the tightening parameters that guarantees 
recursive feasibility, suboptimality, and closed-loop stability. 
Section |VT] presents numerical results using an academic 
example. Finally, Section VII concludes the paper. 




Notation. For u € K”, ||rt|| = ■\/{u, u) is the Euclidean 
norm and [u]^ is the projection onto non-negative orthant 
K”. Given a matrix A, [A]i denotes the i-th row of A and 
j the entry {i,j) of A. Furthermore, 1„ is the vector of 
ones in M" and /„ the identity matrix in In addition, 

eigmax(^) eigjai„(A) denote the largest and the smallest 
(modulus) eigenvalues of the matrix A, respectively. P G 
S"g denotes that P G is positive definite. 


II. Problem Formulation 


Consider the discrete-time linear system described below; 

x{t + 1) = Ax{t) + Bu{t) Vf > 0, (1) 

where x{t) G X C M" denotes the state of the system 
and u{t) G U C K"* denotes the control input. The sets 
X and U are simple proper convex sets (i.e., convex sets 
that contain the origin in their interior). Our goal is to steer 
x{t) to the origin and satisfy the plant constraints. We use 
MPC to achieve these objectives. In this respect, consider 
the following finite-time optimal control problem; 

1 N-l 

V*(a;init) =min - {x'^Qx^+v^Ru^) Xx^PnXn (2a) 

X.U Z 

sX.: x^J^l = Ax^ + f = 0, ...,A^—1 (2b) 

Cx^ + Du^ + g < 0, t = 0,... ,N — 1 (2c) 

To — Tinit (2d) 

Tat G Ajv- (2e) 


where xt and ut are more compact notations for xit) and 
u{t), respectively. For t — — 1 {N denotes the 

prediction horizon), the states and the control inputs are 
constrained in the polyhedral set described by ( [2^ , where 
C G D G g G Note that ( |2cl ) can include 

constraints on the state only, i.e., xt G X, and/or constraints 
on the control inputs only, i.e., Ut G U. In Q G S^o 
and R G S™o- Our problem formulation considers also a 
terminal cost xJ^P^x^ associated with a terminal polyhedral 
setA-AT ;= {x GR^\FNX<fN,FN GRP^^'^Jn GRp^}. 

Through the remaining of the paper, we assume; 


Assumption 1. The pair (A, B) is stabilizable. 

Assumption 2. Suppose Assumption holds. Given the 
gain Kf G obtained by the infinite-horizon linear 

quadratic regulator (IH-LQR)—characterized by the matri¬ 
ces A, B, Q, and R—the following holds: 


Vx e Xn 


X G X, KfU G U, and 
[A BKf)x G p-Xn, 0 < p < 1. 


In addition, the terminal penalty Pn G S"g in the stage 
cost ( |2a| ) is defined by the solution of the algebraic Riccati 
equation associated with the IH-LQR. 


In general, the MPC controller solves the optimization 
problem every time new measurements are available from 
the plant and returns an optimal sequence of states and 
control inputs that minimizes the cost (|^. Fet the optimal 


sequence be defined as follows; 

{x, u} ;= {xg, ...,x*j^,Uq,..., (3) 

Only the first element of u is implemented in closed-loop, 
i.e., the control law obtained using the MPC controller is 
given by; 

KuPcixinit) = Uq, (4) 

and the closed-loop system is described by 

x{t -I- 1) = Ax(/) -I- Bnupcixinii). (5) 


A. Parallelization 


We aim to solve Problem (|^ in parallel. Hence, we exploit 
a similar approach as the one proposed in [3]. Specifically, 
as in [3], Problem (|^ is decomposed along the length of the 
prediction horizon N into N -\- 1 independent subproblems 
to be solved by W -f 1 parallel workers H* (f = 0,..., N). 
Each Ht is allowed to communicate with its neighbours 
Ilt-i and Ht+i at predefined synchronization points. The 
decomposition is possible thanks to the introduction of N 
auxiliary variables Zt {t — 1,..., N) used to break the 
dynamic coupling that arises from ( |2b] l. These Zt can be seen 
as the global variables of the algorithm. In particular, each 
Zt stores the local predicted state Xt+i of each subproblem 
and exchanges this stored information to guarantee consensus 
between neighboring subproblems, i.e., to ensure that the 
predicted state of the (/)-th subproblem, namely is 

equal to the current state of the (/-|-l)-st subproblem, namely 
. Specifically, by introducing the consensus constraints 
zt+i = x[% = defining y* ;= [xf^^ 

Hi := [In 0], H 2 := [A B], p > 0 Problem ([^ becomes: 


N 

min y^Vt(yt,Zt) (6a) 

y,z ' 

t^O 

s.t.: GtVt + gt <0, t = 0,...,N-l (fib) 

Hiyo = Tinit, (6c) 

FIiyN £ (fid) 

zt+i=H 2 yt, t = 0,...,N-l, (fie) 

zt+i = iTiyt+i, / = 0,..., A - 1, (6f) 


where, defining ^t ■= 

• Vo(Co) := where 

Qq ■- p TT T 

“ 2^2 P^n 

and Hq := diag{Q, R}. 


Vt(6) := 

Qtiu where, for f = 1, 

...,N 

-1, 


Ho + p^{H^H2 + HjH^) 

PH, 

PH2 

Qt '■= 

-PH, 

P^n 

0 


L 

0 

pin 

VAr(Ctv) 

= where 




-fAl pin 











Parallelization ^ Relaxation ^ Definition tightening ^ Set solver 

equality constraints parameters accuracy 


Fig. 1: Notation and terminology. 


and Hn := diagjPAf, O^xm}- 

Furthermore, Gt and gt vary for each subproblem as 
follows: 

• Gt := [G D\ and gt '■= g, t = 0,..., N — 1. 

• Gn '■= Op„xm] and gM '■= —fN- 

Remark 1. Note that we introduced a quadratic penalty in 
the cost of the form p/2{\\Hiyt - + \\H 2 yt - -Zt+iP), 

according to the ADMM strategy [4]. This penalty has 
no impact on the cost of the original problem (|^, if the 
consensus constraints are satisfied. 

In the following, we introduce the subproblems that derive 
from Problem Let vt+i (f = 0,..., iV — 1) and wt (t = 
be the Lagrange multipliers associated with the 
equality constraints ( |6^ and ( |^ , respectively. Then, let the 
augmented Lagrangian with respect to the multipliers vt+i 
and Wt be defined as follows: 

£vt+uwt-='Vt{^t)+p[vt+i{H2yt-Zt+i)+wJ'{Hiyt-zt)]. 

Hence, we obtain the following + 1 independent subpro¬ 
blems, called original subproblems, associated with the A^+1 
workers Hj (f = 0,..., N): 

min s.t.. Gtyt T^ 0. (7) 

Vt,zt 

B. Overview of our proposed approach and terminology 

Figure summarizes the main steps that lead to a sub- 
optimal solution of the aforementioned problem and intro¬ 
duces the keywords used in the remaining of the paper. 
Consider the MPC problem (|^ and refer to this problem 
as the original MPC problem. The first step {parallelization, 
detailed in Section |II-A| l is to rewrite the original problem 
in -f 1 independent subproblems (the original subpro¬ 
blems). The aim is to use a Dual Fast Gradient (DFG) 
method to solve these subproblems in order to certify— 
in terms of suboptimality and recursive primal feasibility, 
along with closed-loop stability—the MPC solution. The use 
of the inexact solver will eventually cause a violation of 
the consensus constraints (|6^-(|6f|l introduced to define the 
subproblems Q. Hence, the second step (relaxation equality 
constraints detailed in Section |III-A| l is introduced to relax 
the consensus constraints by a quantitjQ , preventing the 
occurrence of consensus-constraint violations. We refer to 
these subproblems as the equality relaxed (ER) subproblems. 


The set of inequality constraints of the ER subproblems 
includes the inequality constraints of the original subpro¬ 
blems 0 and the inequality constraints due to the relax¬ 
ation of the consensus constraints (|^-(|^. The third step 
(definition tightening parameters detailed in Sections III- 
m|IV[ and|y]l is required to address the following remaining 
issues. Eirst, the solution of each ER subproblem computed 
by the dual fast gradient method might violate the inequality 
constraints, due to the termination of the solver after a 
finite number of iterations. Hence, the constraints of the ER 
subproblems must be tightened by a quantity e* proportional 
to the desired level of suboptimality pt chosen by the 
algorithm. Second, due to the relaxation of the consensus 
constraints, the consolidated prediction, i.e., the predicted 
evolution of the state computed (a posteriori) using the 
control sequence obtained by the independent subproblems, 
might deviate from the predicted local solution computed 
by the independent subproblems and, eventually, violate the 
inequality constraints of the original problem. Hence, an 
additional tightening (dependent on e^t) must be introduced 
on the subset of inequality constraints of the ER subproblems 
that corresponds to the original inequality constraints. The 
proposed algorithm addresses the aforementioned issues by 
exploiting the inequality tightened (IT) subproblems. The IT 
subproblems differ from the ER subproblems in the definition 
of the feasibility region, which is tightened by a quantity 
7t(etj Ezt) that depends on both et and The last step (set 
solver accuracy) selects a suboptimality level pt for each 
subproblem, that guarantees a primal feasible and suboptimal 
solution for the original MPC problem within a fixed number 
of iterations. 


III. Subproblem reformulation 
In the following, we introduce the ER subproblems and 
the proposed algorithm to solve them using At -f 1 parallel 
workers. Eurthermore, we introduce an initial formulation of 
the IT subproblems and derive conditions on the choice of 
the relaxation and tightening parameters to guarantee primal 
feasible and suboptimal solutions for of each subproblem. 

A. Equality constraint relaxation 

Our goal is to obtain a solution for Problem ^ by solving 
the independent subproblems 0 in parallel using inexact 
solvers, such as the Nesterov’s DEG [10]. In order to use 


our proposed solver (introduced in Section III-B 1 , which 
relies on first-order methods, we introduce a reformulation 
of Problem 0 to take into account that the constraints ( |6^ 
and ( |6f| l cannot be satisfied at the equality due to the iterative 
nature of the proposed solver and its asymptotic convergence 
properties. In particular, introducing the relaxation parame¬ 
ters Czt ,Czt+i > 0. for each subproblem (t = 0,..., N — 1), 
the former equality constraints (|6^-(|6^ are replaced by the 
following inequality constraints: 

4 *^ - Zt\<e 


Zt\ f: C-z^ In 'tT' |X/ 


czj In) 


\H2yt ■Zi+ll^Czt+i In 


--Zt-i-il < e 


Zt+1 


(8a) 

(8b) 


'Note that the subscript t indicates that Cz, varies along the prediction Thus, for each subproblem, we can realistically consider a 
horizon. This also holds for the later-defined et, 7 t, and rjt. feasible region defined by the following constraints: 



















.Gt 
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0] 
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9t 

<0, 

(9a) 
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-In 

0] 

6- 
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<0, 

(9b) 

.-Hi 

In 

0] 

6- 

' Czt In 

<0, 

(9c) 

.H2 

0 - 

-In] 

6- 

■ Czt+i Ir 

^<0, 

(9d) 

-H2 

0 

In] 

6- 

Czt+i In 

<0, 

(9e) 


or, in a more compact notation: 

%6+ff6<0. (10) 

where G X(n+m)+2n ._ _|_ 

In the remaining of the paper, we consider the following 
equality relaxed (ER) subproblems: 

V:=minVt(6) s.t.: G^,6 + <?6<0, f = 0,...,iV. (11) 

it 

Hence, let fit ■= [Xf Wt^ f G 

be the Lagrange multiplier associated with the new set 
of inequality constraints defined by where At, wf' 

, v^i, and are the multipliers associated with the 
constraints @, Q, @, and ( |9^ , respectively. 

Then, to handle the complicated constraints ( [TOl i, define, for 
each subproblem, the dual function dt{^t,lJ-t) as follows: 

dtiit.l^t) = min (12) 

where (^t) •— +Mfdiag{/pj,p/4„}(Gjj^t +g^^). 

We refer to ( [T^ as the inner subproblem. Hence, we aim to 
solve the following dual subproblems (outer subproblem) in 
parallel to obtain a solution for Problem (|^: 

df = max dt{^t,lJ't), t = 0,...,N. (13) 

Remark 2. The size of the ER subproblems remains un¬ 
affected if N increases. Intuitively, this modularity is an 
additional feature that can be exploited to preserve some 
favorable numerical properties of the problem (e.g., condi¬ 
tioning, Lipschitz constant, etc.) even when the algorithm is 
running in serialized mode. 

B. Parallel dual fast gradient method 

This section introduces Algorithm that we use to solve 
the MPC problem (|^ exploiting N + 1 parallel workers H^. 
Note that, at this stage, we cannot yet ensure that the com¬ 
puted solution is feasible and suboptimal for Problem Q. 

Algorithm [T] relies on Nesterov’s DEG in which the inner 
problem is solved in an ADMM fashion, as explained below. 
Specifically, as Eigure depicts, at each iteration of the 
algorithrrR Ht computes a minimizer for (steps 

1-4), i.e., the algorithm returns a solution for each inner 
subproblem ( [T^ . In particular, our algorithm, in compliance 
with the ADMM strategy, first minimizes with respect 
to yt in parallel for each subproblem (step 1). Then, using 
the information received by Ht+i, i.e., the updated value 
of Vt+i (synchronization step 2), our algorithm computes— 
in parallel for each subproblem—the value of the global 

^Note that = a;o and to initialize Algorithm 1 


Algorithm 1 Parallel Dual Fast Gradient Method. 

Given , fif, Qt, Gt, gt, Fpt, and kt for each lit (t = 0,..., N) 
while k < kt do 

1. fit computes =argminj^^ C t^tdt , kf) ■ 

1. lit receives Vtii from Ilf+i, send to Ilt-i. 

,fc+l 


3. lit updates ^t ^ according to (ID- 

4. lit receives from Ilt+i, sends 

5. fit computes 


= 


gf + L ^ 


Mt ^ttt 


dtiif+f fit) 


6. Define: a hpt '-—Liit (k+3) 

7. fit computes: 




s=0 


dtii!,fit) 


end while 


Ht Ht+i 




variable zt according to the following rule (step 3): 

4+' = I (ddiyt^^+ff2yt^i+v^+w+-Vt--wf). (14) 


Note that this strategy allows to handle the coupling intro¬ 
duced by the 2-notm in the cost function of ([n). Then, (syn¬ 
chronization step 4) Ht receives (sends) the updated value 
of 4+1^ (4^^) from Ilt+i (to Ht-i), respectively. Einally, 
the worker Hj computes the new values of the multipliers 
p,t~^^ (steps 5-7). We compute offline (for each subproblem) 
the Lipschitz constant associated with Wf^^dti^t, dt) to 


perform the multipliers’ update: 


. = diag 


IIGqII" t \\pdmg{H2,-i}\\l j 

“gmin(Qo)'‘f'“’ eig„,i„(Qo) 2n 


. = diag 


l|Gt 


-Jt 


||pdiag{7/i,-J}||2 


eigmin(Qt) eigmin(Qt) 


||pdiag{7f2,-f}||2 

eigmin(Qt) 


, (f = l,...,V-l). 


LfiN = diag 


IIGjvlls T ||pdmg{gi,-J}||^ y 

“SminlQjv) eig„,.„(Qjv) 2' 


Note that our update rule is different from the one proposed 
in [6], where ADMM is used in combination with Nesterov’s 
fast gradient methods. At each iteration, the algorithm pro¬ 
posed in [6], first computes the exact minimizer yt and then 
updates vt and wt- Our algorithm does not wait until the 
DEG returns a minimizer yt to update the multipliers vt 







































and wt, but starts updating their values along with the DFG 
iterations, encouraging the information exchange between 
neighboring subproblems. This algorithm is also different 
from the one proposed in [3]. In particular, the workers 
exchange the necessary pieces of information before the 
update of the Lagrange multipliers and none of the dual 
variables is exchanged between the neighboring workers, as 
Figure highlights. Furthermore, the information exchange 
between neighboring workers is unidirectional, i.e., Ilt+i 
sends the updated information to flj, but flj does not send 
any updated information to Ilt+i. 

Using an argument similar to the one of Theorem 1 in [8], 
we can compute the primal feasibility violation and the level 
of suboptimality of the solution of each ER subproblem 
returned by Algorithm [T] 

Theorem 1. ([8]) Let Vt(^t) be strongly convex, the se¬ 
quences fl^, fj,^) be generated by Algorithm and 

:= X]s=o (/c+i^)'(a[+ 2 ) T^hen, an estimate on the primal 
feasibility violation for the ER subproblem ([n]i is given by 
the following: 




(15) 


where Rt := \\Pt\\ ■ Moreover, an estimate on primal 
suboptimality is given by the following: 


0<^:-Yt{it)<Rtrif (16) 


Algorithm [T] terminates after a fixed number of iterations 
that depends on rjt and Rt [8]: 


h 


\JsRtr]t 


(17) 


C. Tightening of the original inequality constraints 

In order to guarantee the primal feasibility of each subpro¬ 
blem using Algorithm we introduce + 1 auxiliary sub¬ 
problems, namely the inequality tightened (IT) subproblems, 
which differ from the ER subproblems in the definition 
of the feasible region. In particular, each IT subproblem can 
be defined as follows: 


—rnins.t.: + 56 + ^t'^pt+in 5; 0, (18) 

St 

where ct > 0 is the tightening parameter, which depends on 
the suboptimality level rjt that the proposed algorithm can 
reach within kt iterations ( |T7] i. According to [8], solving ( [TS] ) 
using Algorithm n ensures, with a proper choice of ct, that 
the solution of ( fl^ is primal feasible and suboptimal for 
subproblem O- 

To define an et similar to the one introduced in [8], we 
must compute an upper bound for the optimal Lagrange mul¬ 
tiplier, namely Pt et’ associated with the IT subproblem ( [TSl l. 
We use an argument similar to the one of Lemma 1 in [12]. 
In particular, we compute the aforementioned upper bound 
for pt et according to the following lemma. 

Lemma 1. Assume that there exists a Slater vector jjt G 
such that Gtpt + 5 t < 0. Then, there exists et > 0, 
et < mmj=i^,„^pJ-(Gti/t ezt,ezt+i > £t, such that 
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Fig. 3: Local and consolidated predictions. The different colors highlight 
the different subproblems (e.g., the dark blue color refers to the subproblem 
handled by worker IIo, the blue color refers to subproblem handled by 
worker IIi, etc.). 


the upper bound for pi ^ is given by 


IIm 


t,et I 


< 2i?d, := 2 


- dtjpt) 
min {[Ftjj}’ 

j=l,...,pt+2n'-'^ 


(19) 


where Tt:=[[-{Gtyt+gtV-et ljj[2p(e^,-et) l^][2p(e^,^,- 
et) 1^]]^ G and dt{pt) is the dual function for the 

original subproblem ( [T3] l evaluated at pt G 

Proof: See Appendix [I] ■ 

Remark 3. Lemma does not only provide an upper bound 
for ll^t gj|, but it also provides guidelines to select the values 
of e^t and e^t^t as a function of min^^i -f 

gt)]}, which only depends on the primal variable yt- An 
alternative way to determine the relaxation parameters is to 
include and ezt+t in the set of decision variables and 
penalize them in the cost function as it is usually done 
to handle soft constraints. This will, however, increase the 
number of decision variables in the problem formulation and 
it will have an impact on the original cost. 


IV. Tightening improvement to guarantee primal 
FEASIBLE consolidated PREDICTIONS 


The previous section showed how to choose the tightening 
parameter et of each IT subproblem to ensure that the f-th 
local solution, i.e., the solution computed by the f-th IT sub¬ 
problem ( [T8] |, is primal feasible for the f-th ER subproblem. 
This section provides guidelines to improve the choice of the 
tightening parameter of each IT subproblem •El in order to 
guarantee the primal feasibility of the consolidated solution, 
i.e., the predictions obtained, starting from the initial state 
xq, using the control sequence 
— {p(o) 








( 20 ) 


where the elements of Ug are computed by the independent 
IT subproblems Eigure highlights the difference 

between the local and the consolidated prediction. In par¬ 
ticular, when a new measurement is available from the 
plant, the subsystems ( |T8] l compute in parallel (xo,Uq^^^), 
■ ■ respectively. Ac¬ 

cording to the results of the previous section, the pair 
(xt*^^, ) is primal feasible for the t-th subproblem o, 

thanks to the introduction of the IT subproblems. Never¬ 
theless, due to the relaxation introduced on the equality 







































constraints ( ^aj i -fbll, there is a bounded mismatch between 
and ^ (t = 0,..., N — 1). Hence, starting from the 
initial state xq, when the control sequence Ug is applied to 
compute the consolidated state prediction 

:= {xo,Xi^ei ■ ■ ■ ,XN,eiv}^ ( 21 ) 


the feasibility of Xj is no longer guaranteed. Note, however, 
that Ue G U := Ui X ... xUm, i-e., is feasible. Hence, 
no additional tightening is needed on the input constraints. 

In the following. Section IV-A defines an upper bound 
on the maximal feasibility violation of Xg. This feasibility 
violation is a consequence of the local relaxations of the 
equality constraints. Then, Section [TV-B introduces sufficient 
conditions to ensure the primal feasibility of the consoli¬ 
dated prediction and provides guidelines for the choice of 
the tightening parameters for each IT subproblem. 


A. Upper bound on the maximal feasibility violation of Xj 

Let Uj and x^ be defined by ( |20l i and ( 0 , respectively. 
Moreover, from ( | 8 al i and ( | 8 b| ), the following holds: 

\x[%!_l-x[%\<2e.,. (22) 

Our goal is to characterize how far the consolidated predicted 
state is from the local predicted state. 

Lemma 2. Let the f-step-ahead consolidated prediction xt^et 
be defined by © and assume that ( |2^ holds. Then, there 
exists at € K, > 0 , such that the mismatch between Xt,tt 
and the state of the f-th subproblem x^fl^ is bounded, as 
follows: 

\xt,et - x^tlt\ < (23) 

Proof: See Appendix [n| ■ 

Remark 4. According to Lemma a possible choice of at 
is the following: ^ 

at-.= 2^\\A^\e,,_^. (24) 

j=o 

B. Tightening parameter selection 

According to Lemma 2 xt^et differs from by a 

quantity bounded by at- Inus, it might violate the con¬ 
straints of the f-th subproblem (|7]) by as much as at, in 
the worst-case scenario. In particular, we must ensure that 
CtXt,et + DtUt^et + gt < 0. Using the computed upper 
bound ( |23| l, the following holds: 

^tx[ -f Dtu\. gj + Pt + jUtlat In +et Ipt < 0 

tm 

CltXt,et + + flt + IC'tlCTt In +£* Ipt ^ Oj 

where jCt] indicates the absolute value of Ct- Recall that 
these mismatches are caused by the use of inexact solvers 
and that at depends on In the following, we provide 
guidelines to improve the choice of et for each subproblem. 
Furthermore, we provide a modified upper bound for the 
optimal Lagrange multiplier associated with the tightened 
subproblems ( [TF| , which considers the additional tightening 
introduced by at. 


Lemma 3. Consider the following IT subproblems: 


V 


* 

It 


= inin Vt(^t) s.t.: +56 + 7 t < 0 , 


( 25 ) 


for f = 0 ,..., A^, where 7 *:= [{\Ct\at In +6 e* 
Consider the assumptions of Lemma [T] and the existence of 
at for all f = 1,..., according to Lemma|^ Then, for each 
subproblem, there exist ct > 0 , CztAzt+i > Ct such that the 
upper bound for the optimal Lagrange multiplier associated 
with the IT subproblems (|25ll is described by 


IK^JI < 2nt := 2 . 

Tat := \\—{Gtyt + 5t)^ ~ {\Ct\at In)”^ — 

et)ll][2p{ez,^, - et)ll]]^ 


etljj[ 2 p(e^, - 


Proof: See Appendix [HI] 


Remark 5. The choice of et (f = 0,..., N) is not unique and 
depends on the choice of Czt (f = 1; ■ ■ ■ j N). For example, 
given at in (|24li, a possible choice of (t = 1,..., N) is: 


€zt <min 


Mtv-ti 


min {-(Gtj/t-f pt)j } 

IIAll ’l + 2f max {Er=il[C^tk*l} 


(26) 


Consequently, the tightening parameters are given by: 

et<^minie 2 ,,e;^,+,, rnin {-(Gtjit + 5t)i}G (27) 

for t = 0,... ,N. This choice implies that first we select 
the relaxation parameters and then we adapt the tightening 
parameters on the original inequality constraints based on the 
choice of for all t=l,, N. An alternative is to fix et 
for the inequality constraints and consequently compute Cz^. 
In general, the choice of the parameters strongly depends on 
the system-state matrix A in Q. 


Remark 6 . In the context of this work. Algorithm de¬ 
scribed in the next section, adapts the above derived pa¬ 
rameters at each problem instance. If we consider a fixed 
tightening scheme, such as the one proposed by [7], et and 
ezt can be computed offline (for all the initial states in the 
region of attraction). 


In the following, we show that by using {x.y,u..y}—u..y 
is the control sequence obtained by solving the IT sub¬ 
problems ( |25| l and x-^ is the corresponding consolidated 
prediction—the inequality constraints of the original MFC 
problem (|^ are satisfied. Consequently, the predicted final 
state is in the terminal set of the original problem. If the 
desired level of suboptimality of Algorithm is chosen as: 

Vt := 6 / 2 , (28) 

then, according to Theorem there exists k, 7 t •= 
[vl-it ++ 1 , 7 *]'^ such that + < Vt < et. 

Using similar arguments as in [ 8 ], the following holds for 










+ 9it + 


I I In Ipt 

l4n 


^ lpt+4n ■ 


Hence, for all j = 1,... ,p(, the following holds 


[C'ta;^ + St + ICtjat 1„ +et Ip^Jj 


< et- 


J + 


Consequently, exploiting the upper bound ( |2^ , for all j = 
1,... ,P(, we have; 


+ DtUt^-yt + St + IC'tlcTt In +et lpt]j ^ £t 

which leads to + gt < 0 Vf = 0,..., -/V, 

where Xt^-y^ is the f-step-ahead consolidated prediction com¬ 
puted using the solution to the IT subproblem ( [25l ) with 
tighening parameter jf 

In summary, this section showed that there exists a choice 
of the relaxation and tightening parameters that guarantee a 
feasible consolidated prediction with respect to the original 
problem (|^. 


V. Suboptimality, recursive eeasibility, and 

CLOSED-LOOP STABILITY GUARANTEES 

In the following, we derive bounds for := 

Vt(x,y, u,,,), i.e., the cost obtained using {x,y,u^}, 
with respect to the optimal cost V* of the original problem. 

Theorem 2. Assuming that there exist et (t = 0,. .., N) and 
Czj (t = l,..., N) selected according to Lemma then the 
following holds: N 

V* < < V* + 2 ^ Tlt^tiu (29) 

where ft ■= e* + max { \ [Ct]j,i\}at. 

Proof: See Appendix |IV| ■ 

Theorem established the level of suboptimality of the 
consolidated prediction with respect to the original problem. 
In particular, the sequence {x.y,u..),} is suboptimal for the 
original problem and satisfies the original inequality con¬ 
straints (including those associated with Xn). 

Recall that for the update of TZt, our algorithm requires 
a strictly feasible vector yt for Hence, every time new 
measurements are available from the plant, our algorithm 
must provide a strictly feasible solution (not necessarily 
optimal) for the first pt inequality constraints of each ER 
subproblem. The following lemma provides guidelines to 
compute yt- 

Lemma 4. Let y^ be defined as y^ := = 

[(4 (4.,Then, 

a feasible j/+ at the next problem instance, is given by: 

= [y-t[2:N+i] ((^ + BKf)xN,^j^ff (30) 

Proof: See Appendix [V] ■ 

We want to show that the cost decreases at each problem 
instance. Using a similar argument as in [13], under Assump¬ 
tion on Xm and ensuring that XN,'yN G (thanks to a 


proper choice of the tightening parameters, as the previous 
section showed), we can show that: 


N N 

J2^t{yt)<J2^t{yt,yJ-Vo{yo,^,) Vyo.^o GiVattr, (31) 

where J^attr is the region of attraction. Hence, from ( [2^ and 
iTj, the following holds: 

^ l29l ^ 

E^*(4*) ^ v*(x+) +E/(7t+,7^+) (32a) 

t^o 

N N 

< T.^t{yt) + (32b) 

t—0 t—0 

prj AT N 

^ E - Vo(2/o,7c) + E/(Tt+,7^+) (32c) 

t—0 t—0 


where /(e(+,;= (27^+y^)7+, using 7t+,7^^+ to 
represent the updated values of these parameters according to 
The inequality above shows that the total cost decreases 
at each problem instance if Xn is defined according to As¬ 
sumption and if the A^-step-ahead consolidated prediction 
lies in the terminal set. Asymptotic stability of our controller 
follows if Vo(yo. 7 o) > Yt^o filt Hence, we can 

modify the update of et and e^^ to ensure that is satisfied. 


Remark 7. A possible choice of e^^ {t = to 

fulfill ( |3^ is the following: 


et < min 
— 




(33) 


= Vn 


/ " 

47V7^^+ yff [l + 2t max { E I [^*^41} 


-1 


Consequently, et can be selected according to ( |27| ) to pre¬ 
serve the definition of the upper bound on the optimal 
Lagrange multipliers given in Lemma 

Algorithm summarizes the main steps needed to obtain 
a stabilizing control law when the original MFC problem is 
solved in parallel using inexact solvers. In particular, note 
that, if the measured state is in X^, from Assumption]^ the 
state and the control constraints are automatically satisfied 
without solving the MFC problem in parallel. 


Remark 8. Steps 17-23 are the only nonparallel ones of 
the algorithm (Algorithm [T] is instead fully parallelizable). 
The main reason is in the adaptive nature of the algorithm 
(see also Remark |^. Algorithm adapts et and every 
time new measurements are available from the plant. A fully 
parallel Algorithm is possible using a fixed tightening 
strategy, in which et and can be computed offline. 


VI. Evaluation 

We evaluated Algorithm using the LTI system described 
in [14]. The system (sampled at 7), = 0.5 s) is described by: 

x(t -f 1) = Ax{t) + Bu{t), h{t) = Cx{t) + Du{f), 










Algorithm 2 MPC with adaptive parallel tightening scheme. 

1: Given A,B,X,U,^n,N 
2: Compute offline: Kf, Pf, Fiq, fr<i. 

3: Measure: initial state at time t = 0. 

4: for t = 0 to A do 

5: Compute offline: Gjj, , Qt, Wt, ct- 

6 : Compute: initial strictly feasible vector yt- 

7: Compute: initial tightening according to Lemma 

8: end for 

9: for t = 0 to oo do 

10: Measure: initial state Xinit- 

11: if Xinit S Afjv then 

12: Compute: u = KfX\„[x. 

13: else 

14: Compute in parallel (Alg. 1): § 0 , 74 , ■ ■ •, ?iV, 7 jv exploiting 1251 . 

15: Compute: u = Ujg. _ 

16: Update: y <— u"*" according to J30l. 

17: for f = 0 to A - 1 do 

18: Update: according to Lemma 3 

19: end for 

20: for i = 1 to A do 

21: Update: tt <— according to LemmaM 

22: Update: ■yt ■<— 7 ^ according to Lemma^ 

23: end for 

24: end if 

25: Implement u. 

26: end for 


where x{t) G X :=\x{t) G | |a;i(t)| < 4(j = 1, 2),Vt > 0|, 
u{t) G U := {^(t) £ ]R^||Mi(t)| < 1 (i = l,2),Vt > o|, 
h{t) G % \= {h{t) G < 1 (i = l,2),Vt > Oj, 

and the quadruple (A, B, C, D) is given by: 


1.09 0.22 
0.49 0.02 


1.22 0.88 
-0.78 -0.34 


1.34 -0.16 

-3.19 -0.56 


1.60 1.01 
-0.68 0.77 


The weighting matrices Q, R, and Pn in the cost ( |2al l and 
the IH-LQR gain Kf sxe selected according to [14]. We 
implemented our design in MATLAB (to tune the controller 
and test the initial design) and in C (to run a performance 
analysis). In particular, in MATLAB, we used the Parallel 
Computing Toolbox™ to assign the computation of Algo¬ 
rithm [T] to 8 parallel workers, given a prediction horizon 
N = 7. Furthermore, we relied on the MPT3 toolbox [11] 
to compute Xjy and the optimal solution of Problem (|^. 
Finally, we compared our design to [8]. 

We considered the following scenario. The initial state of 
the system is xq = [—0.101 — 3.7]^. The total number of 
complicated constraints ( |2^ for the original problem is 90. 
We used ( |2^ and to initialize and et, respectively. 
To update them, we relied on ( [33| l and ( |27| ). The selected 
Xq caused u and h to saturate (12 active constraints). In this 
scenario, the state enters Ajv in 3 steps. 

Figure 1^ shows the mismatch between the local prediction 
and the consolidated prediction for one problem instance. As 
Figure depicts, the mismatch (for both states) is below the 
predicted upper bound at for all the A^ -f 1 subproblems. 

Table [I] compares the proposed technique to [8]. The 
table reports the upper bound k on the number of iterations 
needed to achieve a suboptimal solution for Problem 
and the level of suboptimality rj. The table lists only the 
first four subproblems, which are the most significant due 
to the presence of active constraints in these subproblems. 
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Fig. 4: Mismatch between local and consolidated predictions. 
TABLE I 

PERFORMANCE ANALYSIS OF ALGORITHM [T] AND [8]. RESULTS 
SHOW THE MEDIAN OF 11 EXPERIMENTS. 


Sample 

Iterations (for subproblem) 



time 

k ([8]) 

fco 

ki 

k2 

ks ... 

0 

11 • lO"* (2.85 ms) 

58931 (1.85 ms) 

25 

10 

8 

I 

12-10® (302.84ms) 

18218(0.57 ms) 

3480 

0 

0 

2 

11-10® (267.56 ms) 

2265 ( 0.07 ms) 

0 

0 

0 


Suboptimality Level (for subproblem) 



y ([8]) 

^0 

^71 

m 

m ■ • ■ 

0 

3.48 

1.15 

1.15 

1.90 

2.25 

1 

0.31 

0.51 

1.03 

1.21 

1.44 

2 

0.15 

0.62 

1.25 

1.47 

1.74 


The method in [8] and our new proposed parallel algorithm 
produce a comparable behavior, thanks to an appropriate se¬ 
lection of the tightening parameters. The parallel algorithm, 
however, is able to achieve similar results to those in [8] 
using a smaller number of iterations. The larger values of 
k for [8] are probably caused by the value of the Lipschitz 
constant and by the problem conditioning, which affect the 
convergence requiring a higher accuracy for the solver. In 
particular, in our proposed framework, the DFG is applied 
to simpler problems characterized, in the worst case scenario, 
by a Lipschitz constant maXi=o,...,A{i/it} = 196 and by a 
condition number max;^ .. Ar{/Ct} = 104. In [8], the DFG 
solves a larger problem characterized by a Lipschitz constant 
Ls = 21994 and by a condition number Ks = 7020. Hence, 
the modularity of our approach has positive implications on 
important properties for the convergence of the solver. 

Table lists the time required by the optimizer to return 
a suboptimal solution for Problem 0. To measure the 
performance, we implemented both algorithms in C on a 
Linux-based OS. We noticed that given the small size of 
the problem, running Algorithm in parallel did not result 
in significant speedups compared to our algorithm running 
in serialized mode [9]. Nevertheless, in both cases, we 
registered a speedup (230x) compared to [8]. The modularity 
of the proposed algorithm is beneficial even for problems of 
small size, such as the one considered in this section for the 
comparison with [8]. We expect the benefits to be even more 
pronounced when considering problems of larger size. 













































VII. Conclusions 

We proposed an algorithm tailored to MPC that guarantees 
recursive feasibility and closed-loop stability, when the solu¬ 
tion of the MPC problem is computed using inexact solvers 
in a parallel framework. In particular, our algorithm com¬ 
bines ADMM and DFG methods and relies on an adaptive 
constraint-tightening strategy to certify the MPC law. 

Our numerical analysis shows performance improvements 
compared to state-of-the-art nonparallel techniques [8]. Fur¬ 
thermore, our study shows that, for small-size problems, 
even if the solver is implemented in a serialized mode, 
there is substantial performance improvement with respect 
to the state of the art. We expect further benehts from the 
parallelization when the size of the problem increases. A 
scalability analysis of the proposed algorithm on many-core 
architectures is part of our ongoing work. 


Appendix I 
Proof of Lemma[T] 

This section contains the proof of Lemma presented in 
Section IIII-CI 


Lemma 1. Assume that there exists a Slater vector yt G 
such that Gtyt + gt < 0. Then, there exists 


> 0, e* < inmj=i^,„^p^{-{Gtyt + gt)j}, 

eztjCzt+i > such that the upper bound on nl is given 


by 


IImUII :=2 


min {[FiJj}’ 


where 


r* 


— {Gtyt + gt) — e* Ipt 

2p(ezt - g) In 
— et) In 




and dt{jlt) is the dual function for the original subpro¬ 
blem ( fT^ evaluated at fit G 

Proof: The following inequality holds; 


d{fit)<\t{lt) + gZylAM 

= Vt(ft) -I- Xf'^^{Gtyt + gt + etlpj-f 
+ P'^t,*'!^ { — ddlVt + — Czjln + etln) + 

+ (ddlVt — Zt — Eztln + it'Xn) + 

+ PVf+l,et {-H2yt + Zt+1 - ^zt+1 In + etln) + 

+ pvt+lyH2yt - Zt+i - ^zt+i^n + e^ln) 

^ tiit) + Xl'^^{Gtyt + gt + etlp^)+ ( 34 ) 

+ 2p max{u;-;J,u;+;J}(-e^,l„ -f etln)+ 

+ 2p max{ut-+j^,^, 

where the last inequality takes into account that 

G K!^. Dehne wtl := 
maxjwt and := 

Consequently, using the definition of wf and the 

following holds: 


( 35 ) 


Furthermore, recalling that G M+, G K”, and 

G K", the following holds; 


II [At 




t,et 


]^ll<[AtT, <l<i.,]%+2n. (36) 


Hence, if we compute an upper bound for the vector 
[AtJt we obtain an upper bound for H^teJI- 

Thus, from the inequality (|34li, it follows that: 


<Vt{it)-difLt). 


r Ar.., ] 

T 

— {Gtyt + gt) — ctlpt 



2p(ezt - et)ln 

yt+i,€t. 


2 p(e 2 j_,_j — £t)ln 


(37) 


Notice that choosing ct < min {—{Gtyt + gt)j}, 

£zt I Ezt+i > i-e., according to the assumptions of the 

lemma, the elements of Ft are all greater than zero. 

Thus, using ( |T5| ) and it follows: 


9 ■ 1 {F‘]fj-Il/^UII ^ 

Z 3 = l,...,pt+2n 


< [a:,. 


w 


t,et 


V. 




]r. 


(38) 


Consequently, the upper bound on the optimal Lagrange 
multiplier is given by; 


IImUII<2 


Vt{it) - d{fit) 


i=i 


min {[rt]j}' 

,...,pt+2n 


Appendix II 
Proof of Lemma[2] 

This section contains the proof of Lemma presented in 
Section IIV-AI 

Lemma 2. Let the f-step-ahead consolidated prediction x* 
be defined by ( |2T] l and assume that ( |22l l holds. Then, there 
exists at G K, > 0, such that the mismatch between Xj 
and the state of f-th subproblem is bounded, as follows: 

l^t.et - x\% \ < at. 

Proof: In the following, we omit the dependence from 
e* to simplify the notation. The proof is constructive. For 
f = 0, Xq = Xq°\ For f = 1, xi = Axo + i?Mo = which 
is the 1-step-ahead state computed by the local subproblem 
0, i.e., the subproblem associated to worker Hq. Hence, the 
mismatch between xi and x^^^ is simply given by 

1^1 - S^^^l < 2ezi = «!■ 


For t = 2,..., N, the following holds; 


I- -(2) I I- -(1) 
|X2 — X2 I = |X2 — X2 


;^(1) ;^(2)| 


< |X2 - X^^^l -f |X^^^ - X^^^l 

< + Bui — — Bui \ + 2ez 

< 2(||A||e2j -f £ 22 ) = 0^2, 


\\pU\ < ii[a:j* <i,.j^ii- 















<2(||A 


N-l\ 


+ \\A 

+ ■ • ■ + Ezn) = 

which proves the lemma. 


N-2\ 


Appendix III 
Proof of Lemma[3] 


This section contains the proof of Lemma presented in 
Section |lV] 

Lemma 3. Consider the following IT subproblems: 


V; =niin Vt( 6 ) s.t.: G^t + 9L + 7 t < 0 , 


(39) 


forf = 0 , where 7 *:= [{\Ct\at 1„+et e* llnV. 

Given the assumptions of Lemma [T] and the existence of at 
for all f = 1,..., iV according to Lemma Then, for each 
subproblem, there exist e* > 0 , ezt^^zt+t > ^uch that the 
upper bound on the optimal Lagrange multiplier associated 
with the IT subproblems (|39ll is described by 


:=2 


Vtiit) - dtilit) 

minj=i.,..,p,+2„{[raJj} 


f = 0,... ,7V, 


r 




~{Gtyt + gt) — |C't|at In —ejlpj 

2 p(e^, - et)l„ 

~ Ct)ln 


Proof: This lemma follows from Lemma [T] applied to 
the subproblems ( [25| ). From inequality ([3^ formulated for 
subproblem (|25]l, the following must hold 



T 

—{Gtyt + gt) — Yilat 1 „ —ctipt 



2p(e2t - G)ln 

* 

yt+l,et. 


2-p{Czt+i ~ rt)fn 


< \t{it) - d{fLt). 



Hence, in order to satisfy the inequality above, we can select 
the relaxation parameters and the tightening parameters 
according to the assumption of the lemma for t = 0,..., N, 
i.e., the following must hold: 

1 " 

(i) - min {-{Gtyt+ gt)j}>. max { V |[Ct]j 7 |}at + e* 

2 = 1 
n 

(ii) max { Y] | [Gt]j^i | }«* + Ct > 0 

2=1 

(Hi) > et > 0, 


where yt is a strictly feasible solution for the original t- 
th subproblem. Hence, there exists TZt such that the upper 
bound on the optimal Lagrange multiplier is defined as 
follows: 


\K 


7t I 


Vt( 6 ) - d(fxt) 


Appendix IV 
Proof of Theorem|2] 

This section contains the proof of Theorem [^presented in 
Section |Vl 


Theorem 2. Assuming that there exist et (t = 0, ..., N) and 
Czt (f = 1) ■ • ■) TV) selected according to Lemma then the 
following holds: n 

V* <V^ <V* +2 

i=0 

where ft ■= et + max { Ya=i \ [Gt]j,i\}at- 

Proof: Due to the tightening of the original inequality 
constraints, V.^ > V*, given that, as proved in Section IV-B 


the feasible region of the tightened subproblems is inside the 
one of the original subproblems. 

Recall that the consolidated prediction satisfies the equal¬ 
ity constraints ( |2bl l by construction. Hence, the following 
holds: 

N 

V.y<y^ [max(min Vt(a;t,Ut) + {X,GtXt + DtUt + gt)) 

i=0 

+ (Yd [Y 0 ] 7 t)] 

N n 

<V* + 2^72.t7pt(et + max { |[C't]j 7 |}at). 

where \Ip^ 0 ] 7 t selects the first pt components of the vector 
It- ■ 


Appendix V 
Proof of Lemma|4] 

This section contains the proof of Lemma presented in 
Section lYl 

Lemma 4. Let y^ be defined as y^ := [y^o •■•2 /a7„] = 
[(4 ( 4 . 7 J)]^-’Then, 

y+ at the next problem instance, is given by: 

y ~ {y'r[2'.N+i]{iX)-+ BKf)xj\;^jj^) ] 

Proof: We can use a similar argument as the one of 
Lemma 4.2 in [8], recalling that Gt = C for f = 0,..., TV —1, 
Gt = Fm for t = N, Dt= D for t = 0,...,N -I, Dt = 0 
for t = N, and Xn C X according to Assumption]^ ■ 
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